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We describe a Monte Carlo scheme for the grand canonical simulation study of fluid phase equi- 
libria in highly size-asymmetrical binary mixtures. The method utilizes an expanded ensemble in 

which the insertion and deletion of large particles is accomplished gradually by traversing a series 
of states in which a large particle interacts only partially with the environment of small particles. 
Free energy barriers arising from intcrfacial coexistence states arc surmounted with the aid of mul- 
ticanonical preweighting, the associated weights being determined from the transition matrix. As 
an illustration, we present results for the liquid-vapour coexistence properties of a Lennard-Jones 
binary mixture having a 10 : 1 size ratio. 



I. INTRODUCTION 

Fluid mixtures eomprising two or more particle species 
of disparate sizes are common in soft condensed mat- 
ter [1]. A prime example is a colloidal dispersion to 
which much smaller particles have been added such as 
non-absorbing polymers [2-4] or charged nanoparticles 
[5]. Interest in such systems stems from the fact that 
by judicious choice of the small component, one can po- 
tentially control the equilibrium and dynamical proper- 
ties of the large component, giving rise to a rich assort- 
ment of novel phenomena and material properties [2, 6]. 
Given, however, the wide variety of small particles that 
one might conceivably choose to add. the experimental 
task of characterizing the range of possible behaviour is 
considerable. With this in mind there has been much 
interest in deploying statistical mechanics and computer 
simulation to predict the properties of such mixtures. 

In this paper we shall focus on the problem of obtain- 
ing the equilibrium phase behaviour of models of highly 
size-asymmetrical mixtures. Direct analytical assaults 
on such systems are generally complicated by the dis- 
parity in particle length scales [7]. To make progress, 
a widely practiced simplifying strategy is to try to map 
the true two-component mixture onto a single component 
system comprising solely the colloid particles. These are 
assumed to interact via an effective potential which is 
supposed to represent the net effect of the bare colloid- 
colloid interactions plus the additional interactions medi- 
ated by the small particles. Arguably the most successful 
example of such an approach pertains to particles that 
interact as hard spheres - a situation which can be real- 
ized experimentally to a good approximation in colloid- 
polymer mixtures [8]. Here the effective interaction is 
the celebrated "depletion" potential describing the inter- 
action between two hard sphere colloids immersed in a 
"sea" of small hard spheres [9]. In seminal work. Bob 
Evans and coworkers have contributed much insight into 
this situation by tracing out the degrees of freedom as- 
sociated with the small particles in order to produce an 
explicit expression for the depiction potential parameter- 
ized by the particle size ratio and the volume fraction of 
small particles. This not only provides valuable infor- 



mation on the nature of the colloidal interactions, but 
also serves as a basis for theoretical and simulation in- 
vestigations of the phase behaviour of the effective one 
component system [10-12]. 

Whilst impressive progress has been made in obtaining 
accurate effective one-component potentials, at present 
they are largely limited to underlying interactions of the 
hard sphere form [1]. Moreover, because effective poten- 
tials are usually derived in the limit of low density of 
large particles, there are concerns about their accuracy 
at high densities where many body effects are significant. 
Ideally then, one should like to be able to tackle the full 
two component system and treat arbitrary interactions 
between the particle species. Achieving this analytically 
still seems some way off, making it tempting to appeal 
to computer simulation for help. Unfortunately, simula- 
tions of highly size asymmetric mixtures encounter their 
own problems: the relevant physics is controlled by the 
length scale of the large particles, but attempts to re- 
lax these particles are often frustrated by the presence of 
the small particles. For instance grand canonical Monte 
Carlo simulations - the method of choice for studies of 
fluid phase transitions [13] - suffer an unfeasibly small ac- 
ceptance rate for insertions of large particles. Similarly 
in Molecular Dynamics an impractically small timestep 
is mandated by the need to avoid high energy overlaps 
between large and small particles. 

In this paper we describe a tailored Monte Carlo simu- 
lation scheme that circumvents the principal drawbacks 
of traditional approaches. The essential idea is to treat 
both species grand canonically, but to ease the sampling 
bottleneck for insertions (and deletions) of large parti- 
cles by performing these - not in a single Monte Carlo 
step - but gradually. In practice this is achieved by per- 
mitting the system to traverse (in a stochastic fashion) 
a prescribed set of states (or "stages") that interpolate 
between the limits of a large particle being fully present 
and fully absent from the system. This idea of staged 
insertion has been around for some time, principally in 
the context of chemical potential measurements for dense 
fluids and complex molecules using the Widom formula 
[14 19]. It has been recently revisited in the context of 
optimizing expanded open ensembles by Escobedo [20]. 
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However, to our knowledge it has not been used to cal- 
culate the full phase behaviour of a model asymmetric 
mixture at large ratios of the component sizes. 



into a fully interacting (ic. standard) large particle. In 
this sense the n = M state for a system of Ni large par- 
ticles and the n = state for a system of A^; + 1 large 
particles are equivalent. 



II. METHOD 

In this section wc begin by outlining the statistical 
mechanical basis to the staged insertion method for a 
highly size asymmetric binary mixture. Thereafter we 
discuss implementation issues, taking as an example the 
case of a Lennard-Jones (LJ) mixture. 

A. Statistical mechcinics 



n= 1 2 M-1, 12 M-1, O--- 




Ni= 1 2 

FIG. 1: Schematic showing how each integer value of the large 
particle number Ni is expanded into M stages, each of which 
is indexed by the macrovariable n. 



Consider a binary mixture comprising N particles, Ni 
of which arc 'large' (1) and Ng of which arc 'small' (s), 
all contained in a volume V at temperature T. Particles 
are identified via an index 1 < i < N, and a species label 
'ji = l,s, and we write the internal energy as 

N N 

where 4>~f^,-f- is the pair potential for particles i and j of 
species 7i and 7^, located at position vectors q^, and Qj- 
respectively. 

Let us now augment this system with an additional 
'ghost' (G) large particle having position vector q^. The 
ghost particle is taken to interact normally with other 
large particles, but differently with small particles. To 
deal with this, it is more convenient to associate separate 
indices k and m with the Ni large and Ng small parti- 
cles respectively, and write the interaction of the ghost 
particle as 

Ni 

= ^(t^ll{<lk,(lG) + 4>t\^m,(lG) ■ (2) 
fc=l m=l 

Here (j)]^ describes the interaction between the ghost 
large particle and a small particle. This is modified 
with respect to the standard large-small interaction by 
the dependence on a discrete stochastic macrovariable 
n = ... M — 1. The role of n is to index the stages 
that speciiy the degree of coupling between the ghost 
and the small particles. Fluctuations in n forwards or 
backwards across its range result in the gradual inser- 
tion or deletion of a large particle (Fig 1). To be more 
specific, we let n = correspond to A^; large particles, 
while n = M corresponds to iV; -I- 1. Intermediate values 
of n = 1 . . . M — 1 represent a system of Ni large parti- 
cles plus a ghost particle. Thus transitions n = 1 ^ 
correspond to the deletion of the ghost particle from the 
system, while n = M — 1 — )■ M correspond to it turning 



The internal energy of the augmented system is 

^'dq};, {q}s, (jGi = $ + $G and the associated 'ex- 
panded' [21] canonical ensemble has the partition func- 
tion Z'{Ni,N„ V, T, n) , where 



^' = n n / "'qfe / / dqaeM-P^'], (3) 

fe=l TO=1 

with j3 = l/ksT. In the present work, we shall be con- 
cerned with the measured form of the grand canonical 
(GC) ensemble probability distribution of the fluctuat- 
ing number of large particles, p{Ni\iii, ^s,V,T), where 
Us and m are the chemical potentials of the small and 
large species respectively. This is obtainable from mea- 
surements of the joint distribution p{Ni,n\fii, fis,V,T) 
conducted within the expanded GC ensemble, which is 
defined via a weighted sum of the expanded canonical 
ensemble partition function Z': 

00 

p{Ni,n) ^Y.^' [^(^""' + • 

iV3=0 

Here ~ means up to an arbitrary normalization constant 
and (for brevity) we have omitted combinatorical and 
volume factors. p(Af;|/i;, /i^, V, T) follows from Eq. 4 by 
picking out those macrostates from the expanded ensem- 
ble having n = 0, ie. that correspond to the physical 
states in which no ghost particles are present in the sys- 
tem: 

00 M-1 

v{^i) = ^ E E ^' [^(^'-"^ + ^«'"«)] -^".o (5) 

N,=Q n=0 

where 

00 00 JVf— 1 

^ = E E E ^' e^p [^(^'-"^ + ^«'"«)] -^".0 (6) 
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is the grand partition function. 

In tlic present work wc shall seek to obtain p{Ni, n) at 
state points {iii,fXs,T) for which its form may vary over 
many decades. Variations on such a scale preclude di- 
rect measurements of p{Ni, n) unless special biasing tech- 
niques are deployed to facilitate sampling of the regions 
of intrinsically low probability. To this end we utilize 
multicanonical preweighting [22], specifying a sampling 
distribution 



p{Ni,n\w) p{Ni,n) exp[w{Ni,n)], 



(7) 



where w(Ni, n) represents a set of weights defined on the 
discrete combinations {Ni,n}. As shall be described in 
Sec. II B, these weights are chosen such as to ensure 
approximately uniform sampling on the set. The de- 
sired form of p{Ni) is regained from the measured form 
of p{Ni,n) by first using Eq. 7 to unfold the effects of 
the weights, then picking out those macrostates having 
n = 0. 



B. Implementation for a binsiry Lennard- Jones 
mixture 

In order to illustrate how the above formalism can be 
implemented in practice, we consider the case of a binary 
mixture of Lennard-Jones particles. Pairs of particles 
labelled i and j (having respective species labels 7j and 
7j) interact via the potential 



(8) 



Here £7^7^ is the well depth of the interaction, while cr7.7^ 
sets the range of the interaction based on the additive 
mixing rule cr^^^. = ((T7. + a-yj)/2, where (T7. and a^. 
are the particle diameters. Interactions are truncated at 
Tc = 2.5tT7;7j and we take ai as our unit length scale. 

We shall be concerned with state points in which the 
small particles occupy a relatively small fraction of the 
overall volume and act as a quasi-homogeneous back- 
ground to the large ones. Under these circumstances, 
configurations of small particles can readily be sampled 
using a standard GC algorithm at constant chemical po- 
tential, lis- As is customary (in order to make contact 
with experimental scenarios), we choose Hs to yield a 
prescribed volume fraction, rjl, of small particles in the 
reservoir [23]. Since we seek a quasi- uniform density of 
small particles, we set Sgg = Eig = eu/lQ, which ensures 
that the small particle reservoir fluid lies well above its 
own (liquid-vapour) critical temperature. In the results 
of Sec. Ill we refer to a dimensionless temperature which 
is defined as T* = l/{l3eii). 

For highly size-asymmetric mixtures, a large number 
of small particles are typically found within the cutoff 
radius 2.5ais of each large particle. In order to locate 



efficiently these particles, we partition our cubic simula- 
tion box of volume V = into cubic cells of linear ex- 
tent 2.5ais, and maintain a list of cell occupancies. Simi- 
lar cells structures were employed to identify small-small 
and large-large interactions [24]. 

As described in Sec. II A, a large particle is inserted 
or deleted in stages by modifying its interaction with 
the small particles. Accordingly one must specify in ad- 
vance the form of the ghost particle interaction for each 
stage n. Obvious possible strategies include varying the 
well depth of the interaction, or the range. However, 
we have found that neither of these approaches operates 
very effectively in practice because of the rapid increase 
of the potential for distances less than that of the poten- 
tial minimum. Specifically, particles whose separation is 
such that the interaction energy is small at one value of 
n can incur a very high energy penalty at a neighbouring 
stage. This impacts adversely on the acceptance rate, 
a difficulty which can only be mitigated by employing a 
large total number of stages M. 

A superior strategy circumvents this problem by im- 
posing a minimum on the attractive part of the interac- 
tion potential and a maximum on the repulsive part: 



min((/)js(r), (^m2x) r < ais 



max( 



(9) 



Each stage, n, is thus specified by a pair of parameters. 



and 



^max- The form of > 



(r) for two such stages 



is compared schematically with the full potential 4>is{r) 
in fig. 2. 




FIG. 2: Schematic form of the interaction between the ghost 
particle and a small particle, 0j"'(r) (Eq. 9) for two values 
of n, compared to the full LJ interaction potential between 
large and small particles. 

Once the set of stages has been defined, a Monte Carlo 

scheme for sampling them can be implemented. Given 
a system of Ni large particles and a ghost particle at 
stage n, a proposal is made to perform a transition to 
an adjacent stage, n — >■ n' . This proposal is accepted or 
rejected according to a simple Metropolis criterion 
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Pa 



= min(l,exp[-/3($[?'^ -$g'^) + A«;]) , (10) 



where $g is given by Eq. (2) and Atz; = w{Ni,n) — 
w{Ni, n') is the difference in multicanonical weights in the 
old and new states, the specification of which is discussed 
below. Note that special measures pertain to transitions 
that bring the ghost particle to the end of the range of 
n. Specifically, for a transition n = 1 — >■ 0, the ghost par- 
ticle is completely removed from the system; the reverse 
move entails a new ghost being added at a randomly cho- 
sen location. On the other hand, when a ghost particle 
undergoes a transition n = M — I — > Af , it becomes 
fully coupled to the rest of the system, tpuir) = 0is(r), 
and Ni Ni + 1; the corresponding reverse move en- 
tails nominating a randomly chosen large particle to be- 
come a ghost and setting iV; — ^ iV; — 1. In such circum- 
stances the difference in weights appearing in Eq. 10 is 
Aw = 'w{Ni,n) - w{N[,n'). 

In standard GC simulation, updates that insert 
or remove a particle usually incorporate a factor of 
e^^^^'V/f^Ni + l) (insertion) ov e''^'^' Ni/V (deletion) inpacc 
to yield the correct GC equilibrium distribution. When 
operating in the expanded GC ensemble it is convenient 
(in the interests of obtaining a smooth weight function 
in the expanded space of Ni and n) to set the chemical 
potential /i; = and to ignore the volume and particle 
number factors for the time being. The neglected factors, 
as well as the unfolding of the multicanonical weights (cf. 
Eq. 7) are easily accounted for when extracting the final 
GC distribution from the measured form oi p{Ni,n): 



\ogp{Ni\iii) ~ logp(iV;,n = 0|w =0)-F^W^i (H) 
-w{Ni, n = Q) + Ni\ogV - log(7V(!). 

We turn now to the matter of the choices for the num- 
ber of stages M and the associated values of the stage 

parameters (j)'^}^ and ^max. This is governed by three 
main desiderata : 



(i) The rates for transitions between neighbouring 

stages should be roughly equal (in both directions) 
in order to avoid bottlenecks in the sampling. 

(ii) M should be sufficiently large to ensure a reason- 
ably high transition rate. 

(iii) The number of stages M should not be so large 
that the correlation time of the resulting random 
walk in {Ni,n} is excessive (bearing in mind that 
the time to cover a given number of steps grows like 
the square of the number of steps). 

With regard to (i), as we have chosen to implement 
it, staging solely influences the strength of interaction 
between the ghost large particle and the small particles. 
Hence it does nothing to ameliorate the decrease in ac- 
ceptance rate that accompanies an increase in the large 



particle density - a situation analogous to standard GCE 
simulations of single component fluids. Thus even if the 
effects of the small particles were to be offset equally for 
all Ni, one would still expect the transition rate to fall 
with increasing A^;. In such a situation, one can at best 
aim to avoid bottlenecks in the sampling by ensuring that 
(i) is satisfied locally in {Ni,n}. With regard to (ii) and 
(iii), there is in practice a tradeoff to be realized here 
which (in parallel with satisfying (i)) may necessitate a 
degree of trial and error, although more systematic ap- 
proaches have been considered in the expanded ensemble 
literature [20]. In sec. Ill we consider factors affecting 
the choice for one practical situation. 

As discussed in Sec. II A, the form of p{Ni , n) may span 
many decades of probability and in order to sample it ef- 
fectively, multicanonical preweighting is called for. This 
in turn requires knowledge of a set of weights, w{Ni,n), 
that facilitate the even-handed sampling of regions of 
high and low probability. One choice that ensures this is 
w{Ni, n) w — logp(7V(, n) which results in a sampled dis- 
tribution p{Ni , n) that is approximately flat (cf . Eq. 7) 
[25]. However, since p{Ni,n) is just the distribution that 
we seek, the task of determining the weight function ap- 
pears -at first sight to be circular. Fortunately though, 
the situation is saved by the observation that it is possible 
to build up a suitable estimate of w{Ni,n) from scratch 
via iterative means [28]. The approach we favour for 
doing so is based on the transition matrix Monte Carlo 
(TMMC) method [29-33]. 

TMMC works by monitoring the transitions between 
macrostates and using these to infer their relative proba- 
bility. Once sufficient transition statistics have been col- 
lected, it is possible to construct the entire probability 
distribution. The starting point is the macrostate bal- 
ance condition relating the equilibrium probability of two 
macrostates u and v to the transition rates between them: 



p{u)W{u — )■ w) = p{v)W{v — > u) , 



(12) 



where u and v are taken to represent combinations of 
Ni and n. The equilibrium transition rate, W{u — )• v) 
can be estimated in the course of a simulation by ac- 
cumulating the acceptance probabilities for macrostate 
transitions into a collection matrix, C{u ^ v). For every 
proposed move, u — > u, the unbiased acceptance proba- 
bility, a (calculated from Eq. 10 by assuming Aw = 0) is 
added to the collection matrix thus: 

C(m->w) ^ C{u^v)+a (13) 
C{u -^u) ^ C{u ^u) + {l-a) . (14) 

This happens regardless of whether or not the move is 
accepted. 

The transition rates can be extracted from the collec- 
tion matrix via 



W{u v) 



C{u v) 



(15) 
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where the sum in the denominator on the right hand side 
runs over all possible values of the macrovariable. 

Putting the transition rates into equation (12) yields 
the macrostate probabilities 



p{v) _ W{u — )■ v) 
p{u) W{v u) ' 

from which the multicanonical weights follow as 



(16) 



u'(m) — w{v) 



In 



W{u- 



(17) 



Since the collection matrix is concerned solely with un- 
biased acceptance probabilities, one is free to apply an 
arbitrary bias during the simulation without affecting es- 
timates of equilibrium properties. This feature of TMMC 
can be exploited to provide an automated strategy for ob- 
taining a suitable multicanonical weight function. Start- 
ing with no knowledge of the weight function, one simply 
updates w{Ni,n) periodically via equation (17). This al- 
lows the sampling to gradually extend over the range of 
Ni,n, pushing progressively into regions of ever smaller 
probability [34]. Once the region of interest has been 
adequately sampled, the collection matrix provides an 
estimate of the requisite distribution p{Ni,n) via Eq. 16. 
During the simulation we also sample (in list form [35]) 
the instantaneous values of Ni,Ns,n, together with the 
configurational energy $. This permits extrapolation of 
the results for p{Ni , n) in temperature via standard his- 
togram reweighting techniques [36]. 



III. APPLICATION TO THE LIQUID- VAPOUR 
TRANSITION OF A BINARY LENNARD-JONES 
MIXTURE 

As a test of our method, we have applied it to the 

study of liquid- vapour phase coexistence in a LJ mixture 
having particle size ratio q = dss/'^ii = 0.1 and reservoir 
volume fraction of the small particles r/J = 0.01. The sim- 
ulations were performed for a cubic periodic simulation 
box of side L = 7.5, which for this ryj would correspond 
to Ng 8000 in the absence of large particles. Since the 
coexistence properties of this system are known already 
on the basis of simulation studies using a very different 
approach (previously proposed by one of us [37]), there 
exists a convenient baseline for comparison. 

The choice of the stage parameters 4>^^n ^max was 
guided by the criteria set out in Sec. II B. For small 
values of iV; < 130, only two intermediate stages were 
required (ie. M = 3) to obtain a fairly high transition 
rate. However, in order to maintain a roughly constant 
transition rate across intermediate stages for different Ni , 
it was found necessary to vary (^max linearly as a function 
of Ni between the limits shown in Table I. For Ni > 130 
the overlap of the ghost with large particles becomes the 





Q<Ni< 130 


Ni > 130 


Stage, n 


0min (j^max 


0min 0max 


1 


-0.5 7.5^2.7 





2 


-0.8 20 16 


-0.5 7.5 


3 




-0.8 20 



TABLE I: The stage parameters (^J"in '^^^'^ '/'max (expressed in 
units of eia) as used in the simulations. For < Ni < 130, 
two intermediates stages (n — 1,2) were used (ie. M = 3), 
and ^mix was varied linearly as a function of Ni between the 
limits shown (see text). For Ni > 130, three intermediate 
stages were used (M = 4) with no variation of parameters. 



principal ground for rejecting an insertion, and we chose 
to mitigate this by the introduction of an additional stage 
(assigned to n = 1) with parameters 4'min — ^mlx = 0, 
thus making M = 4. No variation of the other stage pa- 
rameters was deemed necessary in this regime, whose val- 
ues for 4'min ^^'^ ^max are included in Table I. Across the 
entire range of TV; studied, the acceptance rate for transi- 
tions varied from ~ 30% at small densities of large parti- 
cle to ~ 5% at liquid-like densities. The principal source 
of this variation is overlaps between the ghost particle 
and large particles; its magnitude compares favourably 
with that occurring in grand canonical studies of single 
component fluids over the same density range. 

The simulations were initialised at the temperature 
T* = 1.047, close to the known critical temperature of 
the model [37]. At this temperature the TMMC method 
was used to obtain a suitable form for the multicanonical 
weight function and thence an estimate of the histogram 
p{Ni,n) for A^; = [0 : 300]. This histogram was then 
reweighted in fxi such as to satisfy the equal area crite- 
rion [38] for the two peaks in the near-coexistence form 
of p{Ni), thereby yielding an estimate of the coexistence 
value of Subsequently the data was extrapolated to 
the lower temperature T* = 1,0 by means of histogram 
reweighting. The resulting form of p{Ni , n) provided a 
suitable multicanonical weight function for a new run at 
this lower temperature, which was again performed for 
rjl = 0.01 (which necessitated a re-tuning of /Xg). By it- 
erating this process we were able to step along the coex- 
istence curve without the need to ever recalculate a mul- 
ticanonical weight function from scratch. Further details 
of this strategy for mapping liquid-vapour coexistence 
lines are described in rcf. [35]. 

Fig. 3 presents the resulting estimates of the coexis- 
tence forms oi p{pi) with [pi = Ni/V) at various tem- 
peratures. Not surprisingly, the distributions exhibit be- 
haviour which is qualitatively similar to that of a single 
component fluid [13]. An estimate of the corresponding 
liquid- vapour binodal can be extracted from the distribu- 
tions (by averaging the density under each peak) and is 
shown in Fig. 4(a). The results are fully consistent with 
unpublished data (to be presented elsewhere) obtained 
using the quite different simulation method of ref. [37]. 
Also included in Fig. 4(a) is the binodal for the single 
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FIG. 3: (a) Estimates of the coexistence forms of P{pi) 
for r)l = 0.01 obtained using the methods described in 
the text. Data are shown for T* = 1.047 (criticality) , 
1.0, 0.95, 0.90, 0.85. (b) The same data expressed on a log 
scale. 



component LJ fliiid determined in a previous study [13]; 
the comparison reveals that the presence of the small 
particles in the mixture depresses the critical tempera- 
ture significantly. Estimates of the phase boundary in 
Hi — T space are shown in Fig. 4(b). 

We point out that obtaining this phase diagram in a 

reasonable timescalc would not have been feasible with- 
out the staged insertion/deletion approach. Our tests 
show that the wall clock correlation time in the absence 
of staging is too large to be reliably estimated. Never- 
theless, a lower bound on the ratio of correlation times 
with and without staging can be estimated via a compar- 
ison of the transition acceptance ratc;s. For r]l = 0.01, 
the insertion/deletion rate without staging is ~ 10~^ at 
liq\iid-likc densities of the large particles. This very low 
acceptance rate is of course attributable to the high like- 
lihood that a randomly chosen large particle insertion 
results in overlaps with one or more small particles - a 
visual impression of the difficulty is provided by config- 



FIG. 4: (a) Coexistence densities (circles) as determined from 
the peak positions of Fig. 3; dots interpolate between the mea- 
sured coexistence densities, and are determined via histogram 
reweighting. Squares show the binodal for the single compo- 
nent LJ fluid obtained in Ref. [13]. Critical points are marked 
(*). (b) Corresponding coexistence points in the Hi —T plane, 
with additional points (dots) obtained via histogram extrap- 
olation. 



urational snapshots of the coexisting phases as shown in 
Fig. 5. Use of staging increases the transition acceptance 
rate to ^ 10~^ for M = 4 stages. The cost overhead is 
an increases in the (round trip) random walk length in 
Ni by a factor of M, thereby increasing the correlation 
time by a factor Af^ ^ 10. Hence we believe that in the 
present case our method is more efficient than standard 
grand canonical sampling by a net factor of ~ 10^. 

Notwithstanding the impressive scale of this speedup, 
the net computational expenditure incurred by our study 
remained significant. This is primarily due to the large 
number of small particles in the system, even for the rel- 
atively low volume fractions of small particles that we 
considered. To be more quantitative, the task of ob- 
taining the initial multicanonical weight function con- 
sumed about a week of CPU time on a 32-core 3 GHz 
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machine, while data collection for each subsequent coex- 
istence state point also took about a week. 




FIG. 5: Configuration snapshots of the coexisting vapour 
phase (upper panel) and liquid phase (lower panel) at T* = 
0.95. 



IV. CONCLUSIONS 

In summary, we have described a grand canonical 
Monte Carlo simulation scheme for the study of fluid 
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